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We consider the three-dimensional Ising model in a L± x x cuboid geometry with finite 
aspect ratio p = L±/L\\ and periodic boundary conditions along all directions. For this model the 
finite-size scaling functions of the excess free energy and thermodynamic Casimir force are evaluated 
numerically by means of Monte Carlo simulations. The Monte Carlo results compare well with recent 
field theoretical results for the Ising universality class at temperatures above and slightly below the 
bulk critical temperature Tc. Furthermore, the excess free energy and Casimir force scaling functions 
of the two-dimensional Ising model are calculated exactly for arbitrary p and compared to the three- 
dimensional case. We give a general argument that the Casimir force vanishes at the critical point 
for p = 1 and becomes repulsive in periodic systems for p > 1. 



I. INTRODUCTION 

The spatial confinement of a fluctuating and highly 
correlated medium may cause long-range forces. The 
Casimir effect, which was theoretically predicted in 1948 
by the Dutch physicist H. B. G. Casimir [T], is a promi- 
nent example. This quantum effect is caused by the 
vacuum fluctuations of the electromagnetic fleld and has 
been proven experimentally in the late 1990's [21[3] . It be- 
comes manifest in an attractive long-range Casimir force 
acting between two parallel, perfectly conducting plates 
in an electromagnetic vacuum. 

Another example for a fluctuation-induced force be- 
ing of similar nature as the Casimir force in quantum 
electrodynamics can be found in the physics of critical 
phenomena |4j|5]. This thermodynamic Casimir effect is 
caused by the spatial conflnement of thermal fluctuations 
near the critical point of a second order phase transition. 
Experimentally, it has been proven for the flrst time by 
measuring the fllm thickness of superfluid '*He fllms as a 
function of the temperature in the vicinity of the lambda 
transition |5J |7] . Since then, the thermodynamic Casimir 
effect was measured in several different systems including 
binary liquid mixtures [8 10_ and tricritical ■^He-'*He 

For several years the shape of the flnite-size scaling 
function determined by Garcia and Chan [6 from the ex- 
perimental data has not been understood theoretically, 
in particular its deep minimum right below Tc. While 
the value of the Casimir force at criticality as well as 
the decay above could be calculated using field theory 
|12| - fr5] . no quantitative results were available for the scal- 
ing region T < Tc except for mean-field-theoretical ap- 
proaches |16[ [TT] . Analytic results exist only for the non- 
critical region below Tc, where contributions to the ther- 
modynamic Casimir force from Goldstone modes [18H20| 
and from the excitation of capillary waves of the liquid- 
vapor ''He interface [5T1 become dominant. 

This unsatisfactory situation was resolved in Ref. [22J, 
where a method was proposed to calculate the thermody- 
namic Casimir force for 0(n)-symmetrical lattice models 
using Monte Carlo simulations without any approxima- 
tions, in contrast to, e. g., the stress tensor method used 
by Dantchev and Krech [23'', which furthermore was re- 



stricted to periodic systems. The Monte Carlo simula- 
tions were done for the three-dimensional (3D) XY model 
on a simple cubic lattice with film geometry T_l x Ly x Ly 
and open boundary conditions along the _L-direction, as 
this system is known to be in the same universality class 
as the superfluid transition in ^He and thus displays 
the same asymptotic critical behavior. The results were 
found to be in excellent agreement with the experimental 
results by Garcia, Chan and coworkers O [7] and for the 
flrst time provided a theoretical explanation for the char- 
acteristic shape of the flnite-size scaling function and in 
particular its deep minimum below Tc. In the following, 
this method was used to determine Casimir forces in vari- 
ous systems and geometries [241 [5S] , while other methods 
for the evaluation of thermodynamic Casimir forces us- 
ing Monte Carlo simulations have also been presented 



In the present work, this method is used to derive the 
universal flnite-size scaling function of the excess free en- 
ergy and thermodynamic Casimir force as functions of 
the aspect ratio p = L\^IL\\ for the 3D Ising model with 
cuboid geometry and periodic boundary conditions. Here 
p is allowed to take arbitrary values from p — > (fllm 
geometry) to p — >■ cxd (rod geometry), while former in- 
vestigations were either at p = |12f[r7] or limited to 
the case p ^ 1 |22ff5n] . The paper is structured as 
follows: In the remainder of Sec. [T] the basic principles 
and deflnitions concerning the thermodynamic Casimir 
effect are discussed and the Monte Carlo method will be 
revisited. In Sec. [H] our Monte Carlo results are dis- 
cussed and compared to recently published results by 
Dohm [3T] who calculated the Casimir force within a 
minimal renormalization scheme of the 0(ri) model at 
flnite p, covering temperatures below and above Tc, as 
well as to fleld-theoretical results obtained for T > Tc 
in the framework of the renormalization group-improved 
perturbation theory (RG) to two- loop order [HI I15| . In 
Sec. |III[ we present an exact calculation of the excess 
free energy and Casimir force scaling functions for the 
two-dimensional (2D) Ising model with arbitrary aspect 
ratios p. We conclude with a discussion and a summary. 
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A. Basic principles 

When a thermodynamical system in d dimensions such 
as a simple classical fluid or a classical n- vector magnet is 
confined to a region with thickness L± and cross-sectional 
area -^jj^^- its total free energy F becomes explicitly size- 
dependent. Then the reduced free energy per unit volume 



foo{T) + Sf{T,L^,L» 



(1) 

can be decomposed [3^ into a sum of the bulk free energy 
density /oo and a finite-size contribution 6f. As we as- 
sume periodic boundary conditions in all directions, the 
surface terms in _L and || directions as well as edge and 
corner contributions are omitted in ([T|). In this case the 
residual free energy Sf equals the excess free energy /ox, 
and we will use /ox instead of Sf in the following. 

In terms of /ox the reduced thermodynamic Casimir 
force per surface area in _L direction is defined as ITS] 



(2) 



where /3 = 1/kBT, and the derivative is taken at fixed 
L||. We omit the index _L for the Casimir force, as we 
will not consider Casimir forces in parallel directions in 
this work. 

When in the absence of symmetry breaking external 
fields the critical point is approached from higher tem- 
peratures, which corresponds to the liquid-vapor critical 
point in the case of a simple classical fiuid or to the Curie 
point in a classical n- vector magnet, the bulk correlation 
length ^oo(i) grows and diverges as [33] 



(3) 



with correlation length exponent i^, reduced temperature 
t — T/Tc^ca — 1, and non-universal amplitude ^_|_. In this 
work we use £,+ = Q^/+ = 1.000183(2) x 0.506(1) valid 
for the 3D Ising model on a simple cubic lattice |34l 135] . 

According to the theory of finite-size scaling [3H] and 
under the assumption, that long-range interactions and 
other contributions irrelevant in the RG sense are negli- 
gible, as for instance subleading long-range interactions 
[37], the thermodynamic Casimir force in the regime 
i-Lj^||)^oo ^ 0,1 where a is a characteristic microscopic 
length scale such as the lattice constant in the case of a 
lattice model, obeys a finite-size scaling form 

(iFc{T,L^,L\\)^ Ll'^d^ix^.p), (4) 

where the scaling variable x±^ can be defined as 

-.^ti'^^\-A^^\ (5) 



p = h\^IL\^\^ denotes the aspect ratio, and 'd}_ is a finite- 
size scaling function. Note that in this work d always 



denotes the scaling function of the Casimir force in _L di- 
rection, while the index describes the reference direction 
_L or II of length L. 

An analogous finite-size scaling relation holds for the 
excess free energy. 



/ox(T,Li,L||) ^LI'^ei(xj_,p) 



(6) 



and -d is related to according to [3T] 

1 x\^d pd 



■&A.{XA.,P) 



1 



v dx I 



dp 



Q^{x^,p). (7) 



The dimensionless finite-size scaling functions and 
are universal, that is, they only depend on gross proper- 
ties of the system such as the bulk and surface universal- 
ity classes of the phase transition, the system shape and 
boundary conditions, but not on its microscopic details 
[371i38j. 

At the critical point To the thermodynamic Casimir 
force becomes long-ranged and for sufficiently large val- 
ues of the length Lj^ asymptotically decays as 

/3^c(Tc,i±,i|i)^LI''?9i(0,p) 

^Ll'^[(d-l)A^(p)-pAl(p)],(8) 

where Aj^(/9) = 0j^(O,p) is the so-called Casimir ampli- 
tude |4j , being - like the finite-size scaling function j_ - 
an universal quantity. Note that for finite aspect ratios 
p > the Casimir amplitude becomes p-dependent. The 
film geometry is recovered by letting p — 0, and Eq. ([s]) 
simplifies to 



/?J-c(To,Li,oo) - LI''(d- l)Ai(O). 



(9) 



Since the 1990s, such universal quantities have been 
subject of extensive theoretical research. They were stud- 
ied by means of exactly solvable models [131 [231 1371 1351 - 
42 , Monte Carlo simulations [9] flG l [22H291 143] . as well 
as within field-theoretical approaches |12f[T51 13T1 HiHlS] . 



B. Reformulation for arbitrary p 

The formulation of the Casimir force finite-size scaling 
laws in the previous section was done by assuming film 
geometry p <C 1, i.e., having in mind the limit p — )■ 0. 
However, if p > 1 this picture is not appropriate and 
should be replaced by a more general treatment. In the 
following we rewrite the basic scaling laws in terms of the 
system volume V = instead of the film thickness 

L^. The resulting scaling functions can be used in the 
whole regime Q < p < oo. 

Using the substitution V p"^^^ in Eq. ([61) we get 



/ox(T,Li,iii)^y-ie(x,p) 



(10) 



with an universal scaling function Q and the generalized 
scaling variable 



x = t 



(11) 
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while the scaling function Qj^ from Eq. ^ is recovered 
ei(xi,p)=/-ie(a;,p), 



as 



with 

i_ 1^ 

Xj_ = p" X. 

Similarly, the Casimir force obeys 

from which we derive the scaling identity 

1 xd pd 



d{x,p)^- 



e{x,p). 



(12) 
(13) 
(14) 

(15) 



dv dx dp 

Note that this identity is equivalent to but simpler than 
Eq. Q. At criticality we now define the generalized 
Casimir amplitude as 

A(p)-e(0,p) (16) 

and find 

/ex(Te,Li,L||)^y-lA(p), (17) 

/?J-c(Tc,Li,L||) ^ -l^-VA'(p). (18) 
The case p — I deserves special attention: As 
d 



dp 



eix,p) 



= 



(19) 



in cubic geometry (see Appendix [A|) , Eq. (15 1 simplifies 
to 

^{x,l) = -^^Q{x,l) (20) 
ai^ ox 

at p = 1, and gives a remarkably simple connection be- 
tween the Casimir force and the excess internal energy, 
Eq. (30), in the cube shaped system, namely 

pFc{T,L,L)^ ^u,^{T,L,L). (21) 

Obviously, the Casimir force vanishes at the critical point 
if p = 1, i. e.. 



79(0,1) = 0. 



(22) 



For completeness we also give the definitions of the 
scaling functions in terms of L|| . As 

/cx(T,Li,L||) ^L^'*e||(a;||,p), (23) 

we find 

e||(x||,p) = p-ie(x,p) (24) 

with Xji = <(Lj|/^+)-'^/'^. Note that in this representation 
the scaling identity Eq. ^ reads 

pjy 

dp 

and in the limit p — > cx) simplifies to 

(xy , oo) = -e||(a;||,oo), 
leading to the simple relation 

/3J"c(r,CX),L||) /ex(r,(X),L||). 



0||(a:||,p) 



(25) 

(26) 
(27) 
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Figure 1. (Color online) Bulk internal energy density eoo(T) 
obtained from three different methods: low temperature series 
[47| (red dotted line), scaling expansion |48| (black line), and 
high temperature series [49| (blue dashed line) . 



C. Monte Carlo method 

In this work we focus on the three-dimensional 
isotropic nearest neighbor Ising model on a L±^ x L\\ x L\\ 
simple cubic lattice with periodic boundary conditions 
and Hamiltonian 



(28) 



where K = j3J > Q is, the ferromagnetic reduced ex- 
change interaction and = ±1 are one-component spin 
variables at lattice sites i. The Monte Carlo simulations 
were done using the Wolff single cluster algorithm |,50] . 
Measuring the reduced internal energy density 



u{T,L^,L«) 



(29) 



the excess free energy and Casimir force is calculated 
as follows [55]: First we determined the excess internal 
energy 

Uex(7,L_L,L||) = u(T, L_L,L||) - Uoo{T) (30) 

by subtracting the reduced bulk internal energy density 
UooiT). We used three different results to get precise 
estimates for Uoo{T) of the 3d Ising model in the differ- 
ent temperature regimes: For low temperatures K > 1/2 
we used the low temperature series expansion (54th or- 
der) by Bhanot et al. ^\ , while for K < 1/6 the high 
temperature series expansion (46th order) by Arisue and 
Fujiwara |49| was utilized. Finally, in the vicinity of the 
critical point we used the expansion recently obtained by 
Feng and Blote flH] , where we also took the bulk critical 
indices [51. = 0.22165455(3), ly^^ ^ 1.5868(3) and 
uj — 0.821(5). These three estimates of Uao{T) show a 
broad overlap, see also the discussion by Feng and Blote 
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|48| . the resulting non-reduced bulk internal energy den- 
sity eoo(T) = Uoo(T)/^ is depicted in Fig. [l] With the 
identity 



we determined /ex by numerical integration, using the 
fact that Ucx goes exponentially fast to zero above Tc 

To obtain the Casimir force, we first calculated the 
internal Casimir force 



PFi{T,L^,L\\) = ^ ^ ^' 32 

which is defined similar to Eq. ^ , by numerical differen- 
tiation, using thicknesses L'^ = L± ± f in order to get an 
integral effective thickness Lj_. With Eqs. ([4 35) and the 
hyperscaling relation dv — 2 — a with specific heat expo- 
nent a, it is straightforward to show that this quantity 
fulfills the finite-size scaling form 



- m{T, L^,^ ) ^ ^-'^''L^r'^^''d'Ax±,p), (33) 
with an universal finite-size scaling function 



dx± 



(34) 



This quantity turns out to be very useful in under- 
standing the Casimir force scaling function -d^ {xj_ , p) for 
p — 0, as will be shown in the next section. Finally, the 



thermodynamic Casimir force is obtained from Eq. ( 32 1 
by integration, 

l3Tc{T,Lj_,^)^ - f —(iFi(T,L^,L\\), (35) 
Jt 

where again the exponential decay above Tc simplifies the 
numerical integration. 



II. RESULTS 

A. Casimir force in film geometry p — > 

In Fig. [2] we plot the internal Casimir force 
Eq. ( [32| , for small aspect ratios p = 1/8 and p = 1/16. 
In the limit of film geometry p — )■ we observe strong 
finite-size effects below the critical point [22J, which are 
caused by the influence of the phase transition in the 
d— 1-dimensional system. In this section we will anal- 
yse this influence in detail and show that is directly 
connected to the specific heat of the d— 1-dimensional 
system. We will give the derivation for periodic systems 
where no surface terms occur, as these terms will com- 
plicate the analysis |24) . 



From the scaling identity Eq. ([t]) for p — ^ 0, 

1 a;_Li9 



r°° dr 

/cx(r, L_L,^||) = -y — Wcx(T,i_L,i||) (31) we get 



ly dx I 



(36) 



(37) 



i.e., within the scaling region and for p— )-0 the 
Casimir force can alternatively be calculated without L±- 
derivative [24] . For the internal Casimir force scaling 
function 



1 ix^d 

d-i — 

ly V ox± 



de±{x^,o) 

dx± 



we find the asymptotic identity 

1 - t 



1 



t 

V 



with the excess specific heat 

Cox(7', i_L,i||) = c(T, i_L,L||) - Coo(r) 



(38) 



(39) 



(40) 



and c = dTu/dT as usual. For p — > 0, this quantity 
contains both the bulk singularity 



(41) 



with amplitudes A±, as well as the singularity of the 
laterally infinite film with finite thickness L±^ at ic(i_L) = 
Tc{L±) /Tc — 1, which scales as 



c(T, L_L,oo) A^. 



(42) 



Here, a* denotes the specific heat exponent of the d— 1- 
dimensional system, A'^ are amplitudes, and the factor 
(Lj^/^+)*^"^" ^/"^ guarantees the correct scaling behavior 
for Lj_ — >■ cxD by cancellation of terms containing a* . 
However, as Ccx enters Eq. (39 1 with prefactor t only. 



the bulk singularity at t = is suppressed (as q < 1) and 
Cex is dominated by the singularity from Eq. ( 42 1 , at 



tc(i±) 



(43) 



The location of the critical point was re-analysed from the 
data of Kitatani et al. [52 including corrections to scal- 
ing, as well as from the data of Caselle and Hasenbusch 
|53| . giving the value x*j_ = —1.535(10). This improves 
the value x]^ = -1.60(2) found by Vasilyev et al. ^5] . 
Furthermore, the other terms in (39 1 are 0(1) near x*j_, 
which leads us to the conclusion that the specific heat 
singularity of the d— 1-dimensional film is directly visible 
in the scaling function •d'j_(xj_, 0) around x± = x*j_, 



x^lx^-xir'^' +0(1). (44) 
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Figure 2. (Color online) Internal Casimir force scaling func- 
tion '&'j_{x±,p) for small aspect ratios p — 1/16 and p — 1/8. 
The black thin line is the extrapolation p — > 0, showing a 
logarithmic singularity at x*j_ = 1.535(10) (see text). 



Figure 3. (Color online) Casimir force scaling function 
'd±{x±,p) for small aspect ratios p = 1/16 and p = 1/8. 
The solid line is the extrapolation p — >■ calculated from the 
integrated logarithmic singularity in '&'j_{x±,0). The dotted 
line is the RG calculation of Griineberg and Diehl |T5) . 



From this arguments we conclude that the scaling func- 
tion ■d'j_(x±,0) has a singularity at x*j_ dominated by the 
specific heat singularity of the d— 1-dimensional system, 
with critical exponent a*. In our case, a* = and the 
singularity is logarithmic. This asymptotic behavior is 
included in Fig. [2] as solid line. 

In Fig. [3] we show the scaling function of the Casimir 
force for p = 1/8,1/16, together with the RG results 
of Griineberg and Diehl [15]. The solid line is the in- 
tegrated extrapolation discussed above. We used a cor- 
rection factor (1 + giL'J^), with gi — —4(1), to account 
for leading systematic errors from the discrete deriva- 
tive, which are expected to be oc in periodic sys- 
tems. The inset is a magnification of the minimum, from 
the divergence of ■d'j_{x±=x*j_, 0) the slope of •d±{x±, 0) 
at x*j_ diverges logarithmically. We find a critical am- 
plitude i?j_(0,0) = -0.310(6) (see Tab." 
within error bars with the values i9_l(0, 
[29J as well as i?_l(0,0) = -0.3052(20) [SI]. The zero 
at '(?^(a;™™,0) (solid line in Fig. p| gives the minimum 
position xl'"" = -1.13(5), with i9jJa;T'".0) = -0.360(5), 
while the finite p results are &j_{x'f'^ ^ -1.10(5), 1/16) = 
-0.352(5) and i?_L(a;™" = -0.95(5), 1/8) = -0.340(5). 



IT| , which agrees 
,0) = -0.3040(4) 



B. Casimir force for finite p 

If we increase p to finite values, the Casimir force scal- 
ing function -d^ {xj_ , p) first changes its shape around the 
minimum. The results for p = 1/6 (Fig. |4^) already 
deviate distinctly from the thinner systems, the mini- 
mum below Tc is not so deep anymore, with '!9_L(a;™'" — 
-0.77(5), 1/6) = -0.323(5). These values deviate only 



slightly from the results of Vasilyev et al. ^29J , x 'j 
-0.681(1) and i?j^(x™", 1/6) = -0.329(1), which we at- 
tribute to the larger statistical error in Ref. [55]. 

When the aspect ratio is further increased to p = 1/4 
(Fig. |4|3), the curve has two minima below and above 
which are nearly equal in depth. Note that the results 
for p > 1 /4 are compared to the predictions of Dohm |.31| 
and show similar behavior. For p > 1/4 the minimum 
below Tc vanishes, while the one above Tc remains. This 
is shown in Fig. [4];, where we plot the Casimir scaling 
function for p = 1/2. 

The results for the cube shaped system with p = 1 
are shown in Fig. |4|i |55] . The case p = 1 is quite in- 
teresting, as here the Casimir force at a; = vanishes 
(Eq. ( [22| ) and even becomes positive for p > 1, although 
the system has symmetric, i. e., periodic boundary condi- 
tions. However, this sign change of the Casimir force at 
p = 1 does not contradict the predictions of Bachas |56| . 
as he assumed an infinite system in parallel direction, 
i.e., p = 0. The scaling function i?(a;, 1) has negative 
slope i?'(0, 1) = -e'(0, l)/dv at a; = 0. This behavior is 
in perfect agreement with Eq. (21 1, as the excess internal 
energy Uq-x^(Tc, L, L) is negative for our model. Further- 
more, i!){x, 1) has a second zero at a; = —2.25(5) where 
u{T,L,L) = Uao(T) holds. Fig. j4h shows results from 



both the calculation using Eqs. ( 29]|35 1 (open symbols) 
as well as Eq. (21 1 (filled symbols JTwhere the latter have 
a much better statistics, as no numerical differentiation 
and integration is necessary. 

Finally, in Fig.|5]we depict the Casimir scaling function 
for values of p larger than one. Now we are in rod geome- 
try and use the appropriate scaling variable Ln instead of 
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Figure 4. (Color online) Casimir force scaling function i?x(a;±,p) for several aspect ratios p — {1/6, 1/4, 1/2, 1}. The dotted 
line is the result of Griineberg and Diehl 15] for p = 0, while the dashed lines are the predictions of Dohm [31 . For p — 1 
we also show results from scaling relation Eq. ( |20[ ) (filled symbols), which have much better statistics, as they are directly 
calculated from the internal energy. 



L±. Due to this rescaling, the scahng function i?||(a;||,p) 
converges to a finite Hmit (a;|| , oo) which should only 
slightly deviate from curves for p = 8, just as in the 
inverse case p = 1/8 (see Fig. |3|. In this regime the 
Casimir force is always positive, leading to a repulsion of 
the opposite surfaces. Note that for p = 8 we increased 
the thickness difference for the calculation of the deriva- 
tive in Eq. ^ to L'j_ = L± ± 4, as, e. g., L± ^ 256 for 
Lii = 32. 



C. Excess free energy 

The excess free energy is shown in Fig. [6] for p < 1. An 
interesting feature of these curves is the non-vanishing 



limit for xj_ — > — oo, which means that for fixed tem- 
peratures T <Tc and — > oo the total excess free 
energy V fc-m approaches a finite value. This behavior is 
a direct consequence of the broken symmetry in the or- 
dered phase [57]: In this phase, which only exists in the 
thermodynamic limit below Tc, the Ising partition func- 
tion is reduced by a factor of two, as the system cannot 
reach the whole phase space anymore. This leads to the 
term — In 2 in the total excess free energy of a periodic 
Ising system below 



e(-oo,p) 



In 2, 



(45) 



independent of shape and dimensionality. Note that, 
e.g., for the g-state Potts model this argument directly 
generalizes to 8(— oo,p) = — Ing. Using Eq. (12l, we 
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Figure 5. (Color online) Casimir force scaling function 
(a; II , p) for aspect ratios p > 1, now as function of the proper 
scaling variable a^y. The dashed lines are the predictions of 
Dohm [3T] . 



Figure 6. (Color online) Excess free energy scaling function 
Q±{x±,p) for several aspect ratios p. The dashed lines are 
the predictions of Dohm [31] for p = 1/4,1/2,1, while the 
solid lines are the limits for x± — > — oo, Eq. (46 1. 
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III. EXACT RESULTS IN TWO DIMENSIONS 



(46) 



this limit is shown as thin solid lines in Fig. [6] The results 
are compared to the field theoretical predictions of Dohm 
|31j . we find a satisfactory agreement for positive and 
also for slightly negative values of x±. Furthermore, our 
value A(l) — —0.63(1) for the cube is compatible with 
the value -0.657(30) obtained by Mon [58J. 

The generalized Casimir amplitude at criticality, 
A(p) = e(0,p) (Eq. ([16])), is listed in Tab. [l] for several 
values of p and is depicted in Fig. [7] together with the pre- 
dictions of Dohm [21] (dashed line) as well as the asymp- 
totes (dotted lines). The inset shows Aj_(p) (circles) and 
A||(l/yo) (squares), showing good agreement with these 
predictions for 1/4 < p < 3. 



Table I. Monte Carlo results for the Casimir amplitudes A(p), 
Ap(p) and t?^(0,p), with p =_L for p < 1 and p =|| for p > 1. 
Note that the critical Casimir force changes sign at p = 1. 



p 


A(P) 


Am(p) 







— oo 


-0.155(3) 


-0.310(6) 


1/16 


-39.8(8) 


-0.155(3) 


-0.310(6) 


1/8 


-9.9(2) 


-0.155(3) 


-0.310(6) 


1/6 


-5.7(1) 


-0.157(3) 


-0.30(1) 


1/4 


-2.60(5) 


-0.161(3) 


-0.290(5) 


1/2 


-0.89(2) 


-0.223(4) 


-0.22(1) 


1 


-0.63(1) 


-0.63(1) 


0.000(5) 


3/2 


-0.68(3) 


-0.45(2) 


0.17(1) 


2 


-0.78(3) 


-0.39(2) 


0.25(1) 


8 


-2.86(5) 


-0.357(8) 


0.36(1) 


oo 


— oo 


-0.36(1) 


0.36(1) 



The scaling function Q± of the excess free energy in 
d = 2 is calculated exactly based on the work of Ferdi- 
nand and Fisher [5^ (Note that the term S,Si{n)T^ /2 is 
missing in Eq. (3.36) of this work). Our scaling variables 



-J 



-J 




Figure 7. (Color online) Generalized Casimir amplitude 
A(p) = 0(0, p), Eq. of the Ising universality class in 

d = 3 (red circles, see also Tab. and in d = 2 (Eq. ( 56 1 , 
blue solid line). The dashed line is the prediction of Dohm 
|31| . while the dotted lines show the asymptotes. The inset 
depicts Ax(p) (circles) and A||(l/p) (squares). 
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1/32 
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Figure 8. (Color online) Excess free energy scaling function 
of the d — 2 Ising model for several aspect ratios p < 1. The 
scaling functions for p > 1 can be calculated using Eq. (52 I. 



differ from theirs, we use x± = t{L±/^^) and p = L±/L\\ , 
while they used r = x±/2 and ^ = 1/p as temperature 
and aspect-ratio variables. 

We start from the partition function of the L± x i|| 
isotropic Ising model on a torus |60| . 



1 



Z(T,Lj_,Ln) = -(2sinh2/s:)5 



II X 



with 



2 

X (Z+ + + Z+ ±Z^), (47a) 
^ above and — below Tc, the four partial sums 

L±-l 

^s"^ n (e3^ii'^^"+^ ±e5^ii'^^"+*) , (47b) 



and cosh7j = cosh 2if coth 2A' — cos(Z7r/Lj^). 

For the bulk free energy density of the d = 2 Ising 
model we using Mathematica f61' derived a nice closed- 
form expression not present in the literature yet, namely 



/oo = -ln(2cosh2if) + — 4^^3 
lb 



3 3 
2 ' 2 

2,2,2 



1,1 



(48) 



with k — 2 tanh 2K/ cosh 2K and the generalized hyper- 
geometric function 4F3(-) pT] , 

After some algebra, the scaling function Q± for arbi- 
trary x± and p can be written as 



(49a) 




(49b) 



Figure 9. (Color online) Casimir force scaling function 
i}f_i{Xfi,, p) of the d = 2 Ising model for several aspect ratios p. 
Shown is '&±{x±,p) for p < 1 and i?||(a;j|,p) for p > 1. Note 
that i3±{x±, p) = p'^'&w (i|| , p). 



and 



I±{x^_)= f da;ln(l±e-V^i+^*^ 

J —OO 



Note that 



I±{x±) ^ lim -\nP^{rx±,r) 



(49c) 



(50) 



independent of d. As the 2d system is invariant under 
exchange of the directions _L and 1 1 , 



e(a;,p)-e(x,l/p), 
which using Eq. ( |12[ ) gives 

Q±{x±,p)/p = pQ±{xj_/p,l/p), 
we can derive the identities 



(51) 
(52) 

(53a) 

(53b) 
(53c) 



which are a generalization of Jacobi's imaginary trans- 
formations for elliptic iS functions [62) . 

The resulting excess free energy scaling function 
Qj_{xj_,p) is depicted in Fig. [s] showing a similar be- 
havior as in the three-dimensional case. For a;j^ — )• — oo 







p+,{x^/pMp) 


gl+{x±_)/p 


Py^i.x^,p) 




p+{x^/pMp) 


~ g/+(a;i)/p 


Po{x±,p) 


gp/_(a;j_/p) 


P^{x^/p,l/p) 





Eq. ( 49 ) simplifies to 

e_L(-oo,p) 



-pin 2, 



(54) 
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as explained in Sec. |II C| 

From Eq. ( 49 ) we directly obtain values of the scaling 



function at the critical point xj_ — 0, as 
^+(0) = 



TT 

12' 



(55a) 



and 



^v.(o,p) = (T«;g')L, ^o*(o,p) = ^(Ti;g')L (55b) 



with q = e '^^^ and the g-Pochhaninier symbol [5T] 
(a; q)oo, leading to 



6^(0, p) - -pin 



2qi/i2 



4q-l/6 



^ ''(4^2(0, 9)^93(0, g)^4(0,<z))i/3 ^ ^ 

after expressing the g-Pochhammer symbols in terms of 
elliptic I? functions. This result was already given by 
Ferdinand and Fisher [53] (Eq. (3.37)). The resulting 
Casimir amplitude A(p) is shown as blue solid line in 
Fig-H 

From the exact solution Eq. ( [49| we calculated the 
Casimir force scaling function by numerical differentia- 
tion using the scaling relation Eq. ([t]), as an analytic 
derivation would be too lengthy for arbitrary p. The re- 
sults are shown in Fig. |9] for p < 1 we show (x_l , p) , 
while for p > 1 we show z^y (a^y , p). Clearly the Casimir 
force changes sign from negative to positive values with 
increasing aspect ratio p, as in the three-dimensional 
case. 

Finally we give expressions for the limits p — > and 
p — > oo. In film geometry, p — > 0, Eq. (49 1 reduces to the 
simple result 



2 



dtj In 1 



(57) 



yielding the already exactly known Casimir force scaling 
function 1631 



1 



(58) 



In the opposite limit p — > oo we have 

ey(a;y,Oo) = -'!9||(xy,0o) = -/+(xy) (59) 



Table II. Signs of the terms P^{x_i_, p) in Eq. (49a I for differ- 
ent boundary conditions. 



BCx 


BC|i 


p+ p- 






periodic 


periodic 


+ + 


+ 




periodic 


antiperiodic 


+ + 




+ 


antiperiodic 


periodic 


+ 


+ 


+ 


antiperiodic 


antiperiodic 


+ 


+ 


+ 



using Eq. ( 26 ) . For both p — and p — >■ oo we have the 
symmetries 6j^(a;j^,p) = 6_l(— a;_L,p) and i^j. (a;_L , p) = 
t9j^(— .Tj^, p). Note that all scaling predictions from the 
previous sections have been verified in the d — 2 Ising 
case. Finally, we remark that these calculations can be 
easily extended to mixed periodic-antiperiodic boundary 
conditions by mo difyin g the prefactors of the four terms 
Ps'ix^^p) in Eq. (49a I according to Tab.|ll] 



IV. SUMMARY 

In this work we calculated the universal excess free 
energy and Casimir force scaling functions, Q{x, p) and 
i9(a;,p), of the three- and two-dimensional Ising model 
with arbitrary aspect ratio p and periodic boundary con- 
ditions in all directions. In d = 3 we used Monte Carlo 
simulations based on the method by Hucht |21|, while 
in d = 2 we derived an analytic expression, Eq. ( [49| , 
for the excess free energy scaling function Q{x,p). Fur- 
thermore, we derived several new scaling identities for 
the scaling functions: We showed that the Casimir force 
scaling function -d^ {x± , 0) in the film limit has a sin- 
gularity of order {x± — a;^)^~" at the point < 
where the d— 1-dimensional system has a phase transition 



(Eq. (44 1), where a* denotes the specific heat exponent 
of the d— 1-dimensional system. In our case a* ~ and 
the singularity is logarithmic as shown in Figs. |2] and |3] 
At finite values of p > 1/4 our results are compared to 
field-theoretical results of Dohm, and we find good agree- 
ment in the regime a; > — 2 where his theory is expected 
to be valid [3T] . For the cube with p = 1 we observed an- 
other interesting result, here the Casimir force vanishes 
exactly at the critical point, i?(0, 1) = 0. In Appendix [A| 
this property is shown to hold for all systems that are 
invariant under permutation of the directions, and is not 
restricted to periodic systems. The vanishing Casimir 
force could serve as a stability /instability criterion with 
respect to p: If we assume that the system can change 
the lengths at constant volume, we see that the cube 
with p = 1 and periodic boundary conditions is unstable 
under variation of p at a; = 0, as p < 1 tends to p — ?> 
and p > 1 tends to p — >■ oo. Note that this behavior 
would reverse for antiperiodic boundary conditions, then 
the cube would be stable at a; = and the equilibrium 
shape would even be temperature dependent, as the zero 
of 'd{x, p) varies with x, see Fig. [9] For p > 1 the Casimir 
force is positive and converges against the negative excess 
free energy, ??y(a;y,oo) = -ey(a;y,oo), Eq. ([26 1. 

The excess free energy below Tc is /ex ~ — l^^^ln2 
in periodic Ising systems [57] independent of system 
shape (Eq. (45)), leading to a finite p-dependent limit 
of e_L(-oo,p)7Eq. (|46]). 

Finally, the universal scaling function 0_l (a; , p) is cal- 
culated exactly in d = 2, and the results are found to 
be in qualitative agreement with the results for d = 3. 
The most important difference between these two cases 
is the fact that the 2d system has several symmetries 
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under the condition 



Figure 10. (Color online) The d— 1-dimensional plane B of 
constant volume for d = 3, viewed from the normal direc- 
tion (1,1,1). The origin at the center (o) is the cube with 
6 = 0, while the filled symbols are deformed systems as in- 
dicated by the pictures: The black points mark the direc- 
tions zbs), with constant symmetric under permutation V, 
Eq. ( |A3[ ). The blue arrow Sp (p > 1) and the red arrow —Sp 
{p < 1) mark the direction of the shape variation in terms of p 
used in this work. The black curve sketches a line of constant 
Q{x,b). Note that Q{x,Sp) ^ Q{x, —Sp), as the shape and 
thus is not symmetric under the transformation p 1 /p, 
see also Fig. [7] 



not present in the 3d system, i.e. {x±,p) o (2:||,l/p), 
(a;, 0) -n- (— x, 0), and (a;,oo) -f-?- {—x,oo). 



d 



(A2) 



defining the plane B with constant volume . The sym- 
metry under permutation of the d lattice axes implies 



e(x,5) = Q{x,V{h)) 



(A3) 



with permutation operator V . This symmetry holds if the 
boundary conditions in all directions are equal. Without 
loss of generality we now assume d — i, hi — Q and vary 
the shape of the system along directions 2 and 3, i.e., 
1)2 — — ^3, so that Q{x,esi) = Q{x, —e si) with si — 
(0, 1, —1) and real e. Hence 0{x, e si) is an even function 
of e and thus the directional derivative along si at the 
origin vanishes, 



d_ 



Q{x, e si 



0. 



(A4) 



The same argument holds for the symmetric directions 
S2 = (-1,0,1) and S3 = (1,-1,0) (see Fig. [lO|. As the 
d{d~ l)/2 vectors form an (over)complete base in the 
d— 1-dimensional plane B, and all directional deriva tives 
vanish at the origin 6 = 0, we conclude that Eq. (A4| 
holds for all directions s G B. Hence Eq. ( A4 1 also holds 
for the special case Sp = (2/3, —1/3, —1/3) which is the 
direction of the shape variation used in this work (with 
p = e'^), if we set Li = L± and L2 — — iy . From this 



we conclude that Eq. ( 10 1 satisfies 



d_ 
dp 



Q{x,p) 



0. 



(A5) 



Note that these arguments can be extended to weakly 
anisotropic systems, while less is known in the strongly 
anisotropic case |^ [M] . 
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Appendix A: Stationarity of Q{x,p) at p = 1 

The stationarity of the excess free energy scaling func- 
tion Q{xtP) at p — 1 can be derived for isotropic sys- 
tems with arbitrary symmetric boundary conditions and 
in arbitrary dimensions d: We allow arbitrary shape 
changes of /cx(T, Li,...,Ld 
that Eq. (10) now reads 



and write = e f'L, so 



''"L) - L-'^e{x,b), 



(Al) 



Appendix B: Proof of i9(0, 1) = in the large-n limit 

In this appendix we show for the large-n limit |31] 
that the finite-size scaling function of the thermodynamic 
Casimir force vanishes at bulk criticality in the case of 
a cubic system geometry p = 1. To this end we start 
from the scaling function of the singular free energy per 
volume given by Dohm (Eq. (3) in [31j), together with 
the self-consistency equation for the parameter P{xj_, p) 
at x± — 0, 



F(0,p) = -4^gi(P(0,p)2,p), 



(Bl) 



and the functions Gj{P^, p) (Eq. (4) in [31]). Introducing 



the parameter P{p) 
tegration variable z 



p^P{Q, p) and furthermore the in- 
: p^ z in the integral Qo{P^,p), the 
value of the excess free energy scaling function Q± (see 
Eq. (l6|) at bulk criticality can be cast in the form 



e^{o,p)^A^{p) 



,'A{p) 



(B2) 
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upon setting ^ = —2/3 and 5 = 4/3, where A(p) is given 

by 



A(p) 



P{pf , 1 

127r 2 

3/2 



°°dz / 

— exp 



47r2 



(p2/3£)2l . (B3) 



According to Eq. ( 18 1 one has 



i?±(xi=0,p=l) = -A'(l) 



(B4) 



where the derivative of A(/3) with respect to p at p = 1 
becomes 



A'(l) 



P(1)P'(1) 



47r 



P(l)+ 471^1 (P(l)2,l) . (B5) 



Since P(l) = P(0, 1) is the solution to Eq. (Bl I at p 



1, the expression in square brackets vanishes and thus 
A'(l) =0. 
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